R Notebook

FOLLOW-UP – Add preamble to each MakeData script

Goals

Initialize metadata for the GRS analyses Identify individuals who are eligible for analysis based on having known Celiac status and no known first degree relatives in the analysis cohort

Load frequently used packages

library(openxlsx)
library(ggrepel)
library(dplyr)
library(tidyr)
library(data.table)
library(broom)
library(broomExtra)
library(tibble)
library(sjstats)
library(car)
library(lme4)

Loading required package: Matrix
Warning: package ‘Matrix’ was built under R version 4.1.2

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(lmerTest)

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

library(ggplot2)
library(tibble)
library(modelr)
library(tidyverse)
#library(miceadds)
library(ggforce)
require(openxlsx)
library(tidyverse)
library(caret)
library(glmnet)

Loaded glmnet 4.1-3

library(ggplot2)
select <- dplyr::select
filter <- dplyr::filter

Print session info for reproducibility

sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] glmnet_4.1-3                 lmerTest_3.1-3               lme4_1.1-29                  Matrix_1.4-1                 scales_1.2.0                
 [6] IlluminaDataTestFiles_1.30.0 illuminaio_0.34.0            ComplexHeatmap_2.8.0         RColorBrewer_1.1-3           circlize_0.4.14             
[11] lmPerm_2.1.0                 MASS_7.3-56                  gridExtra_2.3                caret_6.0-91                 lattice_0.20-45             
[16] ggforce_0.3.3                forcats_0.5.1                stringr_1.4.0                purrr_0.3.4                  readr_2.1.2                 
[21] tidyverse_1.3.1              modelr_0.1.8                 car_3.0-12                   carData_3.0-5                sjstats_0.18.1              
[26] tibble_3.1.6                 broomExtra_4.3.2             broom_0.7.12                 data.table_1.14.2            tidyr_1.2.0                 
[31] dplyr_1.0.8                  ggrepel_0.9.1                ggplot2_3.3.5                openxlsx_4.2.5               skimr_2.1.3                 

loaded via a namespace (and not attached):
  [1] readxl_1.4.0         backports_1.4.1      plyr_1.8.7           repr_1.1.4           splines_4.1.1        listenv_0.8.0        TH.data_1.1-0       
  [8] digest_0.6.29        foreach_1.5.2        htmltools_0.5.2      fansi_1.0.3          magrittr_2.0.3       cluster_2.1.3        doParallel_1.0.17   
 [15] tzdb_0.3.0           recipes_0.2.0        globals_0.14.0       gower_1.0.0          matrixStats_0.61.0   R.utils_2.11.0       sandwich_3.0-1      
 [22] askpass_1.1          hardhat_0.2.0        colorspace_2.0-3     rvest_1.0.2          haven_2.4.3          xfun_0.30            crayon_1.5.1        
 [29] jsonlite_1.8.0       survival_3.3-1       zoo_1.8-9            iterators_1.0.14     glue_1.6.2           polyclip_1.10-0      gtable_0.3.0        
 [36] ipred_0.9-12         emmeans_1.7.3        GetoptLong_1.0.5     sjmisc_2.8.9         future.apply_1.8.1   shape_1.4.6          BiocGenerics_0.38.0 
 [43] abind_1.4-5          mvtnorm_1.1-3        DBI_1.1.2            Rcpp_1.0.8.3         xtable_1.8-4         performance_0.9.0    clue_0.3-60         
 [50] bit_4.0.4            stats4_4.1.1         lava_1.6.10          prodlim_2019.11.13   datawizard_0.4.0     httr_1.4.2           ellipsis_0.3.2      
 [57] R.methodsS3_1.8.1    pkgconfig_2.0.3      farver_2.1.0         sass_0.4.1           nnet_7.3-17          dbplyr_2.1.1         utf8_1.2.2          
 [64] labeling_0.4.2       tidyselect_1.1.2     rlang_1.0.2          reshape2_1.4.4       effectsize_0.6.0.1   munsell_0.5.0        cellranger_1.1.0    
 [71] tools_4.1.1          cli_3.2.0            generics_0.1.2       sjlabelled_1.1.8     evaluate_0.15        fastmap_1.1.0        yaml_2.3.5          
 [78] bit64_4.0.5          ModelMetrics_1.2.2.2 knitr_1.38           fs_1.5.2             zip_2.2.0            future_1.24.0        nlme_3.1-157        
 [85] R.oo_1.24.0          xml2_1.3.3           compiler_4.1.1       rstudioapi_0.13      png_0.1-7            reprex_2.0.1         tweenr_1.0.2        
 [92] bslib_0.3.1          stringi_1.7.6        parameters_0.17.0    nloptr_2.0.0         vctrs_0.4.0          pillar_1.7.0         lifecycle_1.0.1     
 [99] jquerylib_0.1.4      GlobalOptions_0.1.2  estimability_1.3     insight_0.17.0       R6_2.5.1             IRanges_2.26.0       parallelly_1.30.0   
[106] codetools_0.2-18     boot_1.3-28          assertthat_0.2.1     openssl_2.0.0        rjson_0.2.21         withr_2.5.0          S4Vectors_0.30.2    
[113] multcomp_1.4-18      mgcv_1.8-40          bayestestR_0.11.5    parallel_4.1.1       hms_1.1.1            rpart_4.1.16         timeDate_3043.102   
[120] base64_2.0           coda_0.19-4          class_7.3-20         minqa_1.2.4          rmarkdown_2.13       Cairo_1.5-15         pROC_1.18.0         
[127] numDeriv_2016.8-1.1  lubridate_1.8.0      base64enc_0.1-3     

Define in/out directory

dir <- "/Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data"

Load MEGA to HTP ID key

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

MEGA.IDkey <- fread("MEGA_041822_MEGA2_to_HTP_ID_key_v0.1_JRS.csv")

MEGA.IDkey

Read in latest version of the meta/comorb data

setwd("/Users/shawjes/Dropbox/EspinosaGroup/DATA_MAIN/META/OneDrive_1_2-17-2022")

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/DATA_MAIN/META/OneDrive_1_2-17-2022 inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

meta_comorb.022021 <- fread("HTP_CLEANED_02_2021_v0.5_COMORB_MASTER_RecordID_vs_LabID.Labels.tsv") %>%
  select(FamilyID, RecordID, LabID, Age_at_visit, Sex, Karyotype, everything()) %>%
  mutate(Celiac = as.numeric(`Celiac disease`))

meta_comorb.022021

Put the column names representing identifiers into a vector for reference below

ID_colnames <- c("FamilyID", "RecordID", "MEGA.FID", "MEGA.IID", "MEGA.LabID")

Make subset of meta/comorb dataset specific to Visit 1 for T21s with complete Celiac status

meta_comorb.T21_visit1_Celiac <- meta_comorb.022021 %>%
  filter(Karyotype == "T21" &
           !is.na(LabID) & Event_name == "Visit 1" & # Assuming that Celiac status for HTP was ascertained as of participant visit 1
           !is.na(`Celiac disease`)) %>%
  select(colnames(meta_comorb.022021)[1:18],
         `Celiac disease`, Celiac,
         everything()) %>%
  rename(Age_at_visit1 = Age_at_visit)

meta_comorb.T21_visit1_Celiac

Make meta/comorb dataset specific to visit 1 for T21s with complete Celiac status and complete MEGA genotyping

#colnames(MEGA.IDkey)
#colnames(meta_comorb.022021)

meta.MEGA.T21_visit1_Celiac <- MEGA.IDkey %>%
  filter(HTP_participant == "Yes") %>%
  select(-c(HTP_participant)) %>%
  left_join(meta_comorb.T21_visit1_Celiac %>% mutate(In_T21celiac_metadata = 1),
            by = c("FamilyID", "RecordID"))
# Note: Here I'm joining to the metadata corresponding to visit 1, when Celiac status was determined for the purposes of the study. 
# For this reason I am not including LabID as a join by variable.

meta.MEGA.T21_visit1_Celiac %>%
  group_by(In_plink_fam, In_T21celiac_metadata) %>%
  summarise(N = n())

`summarise()` has grouped output by 'In_plink_fam'. You can override using the `.groups` argument.

print("6 participants have complete genotype data but incomplete Celiac status.")

[1] "6 participants have complete genotype data but incomplete Celiac status."

meta.MEGA.T21_visit1_Celiac01 <- meta.MEGA.T21_visit1_Celiac %>%
  filter(!is.na(In_plink_fam) & !is.na(In_T21celiac_metadata)) %>%
  select(-c(In_plink_fam, In_T21celiac_metadata))

meta.MEGA.T21_visit1_Celiac01

Define custom function for mode

# Create the function.
getmode <- function(v) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}

Check that each individual was only genotyped one time

if( (meta.MEGA.T21_visit1_Celiac01 %>%
     group_by(RecordID) %>%
     summarise(N = n()) %>%
     arrange(desc(N)) %>%
     ungroup() %>%
     summarise(var(N)) == 0) ) {"Great, all RecordIDs have the same number of rows of data."} else {
       print("Oops, RecordIDs have varying number of rows.");
       temp <- meta.MEGA.T21_visit1_Celiac01 %>%
         group_by(RecordID) %>%
         summarise(N_rows = n()) %>%
         arrange(desc(N_rows)) %>%
         ungroup() %>%
         mutate(min.N_rows = min(N_rows),
                max.N_rows = max(N_rows),
                mode.N_rows = getmode(N_rows),
                N_rows_Equal_mode = ifelse(N_rows == mode.N_rows, 1, 0),
                N_rows_NotEqual_mode = ifelse(N_rows != mode.N_rows, 1, 0))
       print(paste("N = ", sum(temp$N_rows_Equal_mode), " RecordIDs with N_rows = ", unique(temp$mode.N_rows), sep = ""))
       print(paste("N = ", sum(temp$N_rows_NotEqual_mode),
                   " RecordIDs with mean N_rows = ",
                   mean((temp %>% filter(N_rows!=mode.N_rows))$N_rows),
                   "\nPrinting problematic RecordIDs now. Problems stored in `problems`",
                   sep = ""))
       #print(paste("N = ", sum(temp$N_Equal_mode[1]), " samples with nrow = ", temp$N, sep = ""))
       
       problems <- temp %>%
         filter(N_rows != mode.N_rows) %>%
         select(RecordID, N_rows, mode.N_rows) %>%
         left_join(meta.MEGA.T21_visit1_Celiac01, by = "RecordID") %>%
         select(FamilyID, RecordID, MEGA.IID, MEGA.LabID, N_rows, mode.N_rows) %>%
         unique()
       
       print(problems)
     }

[1] "Oops, RecordIDs have varying number of rows."
[1] "N = 205 RecordIDs with N_rows = 1"
[1] "N = 5 RecordIDs with mean N_rows = 2\nPrinting problematic RecordIDs now. Problems stored in `problems`"

print("Based on this output, five samples were genotyped twice in different wells.")

[1] "Based on this output, five samples were genotyped twice in different wells."

print("Below we will clean the data to keep only the first sample well for each RecordID.")

[1] "Below we will clean the data to keep only the first sample well for each RecordID."

Since some RecordIDs/LabIDs were genotyped in more than one well, sort MEGA.IID alphabetically within RecordID and then keep the first MEGA.IID

Check if now each RecordID has only one MEGA.IID

print(paste("Rows in `mega_Snp_dosage.clean01` = ", meta.MEGA.T21_visit1_Celiac01 %>% select(RecordID, MEGA.IID) %>% unique() %>% nrow(), sep = ""))

[1] "Rows in `mega_Snp_dosage.clean01` = 215"

print(paste("Rows in `mega_Snp_dosage.clean02` = ", meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>% select(RecordID, MEGA.IID) %>% unique() %>% nrow(), sep = ""))

[1] "Rows in `mega_Snp_dosage.clean02` = 210"

if( (meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
     group_by(RecordID) %>%
     summarise(N = n()) %>%
     arrange(desc(N)) %>%
     ungroup() %>%
     summarise(var(N)) == 0) ) {"Great, now all RecordIDs have the same number of rows of data."} else {
       print("Oops, RecordIDs have varying number of rows.");
       temp <- meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
         group_by(RecordID) %>%
         summarise(N_rows = n()) %>%
         arrange(desc(N_rows)) %>%
         ungroup() %>%
         mutate(min.N_rows = min(N_rows),
                max.N_rows = max(N_rows),
                mode.N_rows = getmode(N_rows),
                N_rows_Equal_mode = ifelse(N_rows == mode.N_rows, 1, 0),
                N_rows_NotEqual_mode = ifelse(N_rows != mode.N_rows, 1, 0))
       print(paste("N = ", sum(temp$N_rows_Equal_mode), " RecordIDs with N_rows = ", unique(temp$mode.N_rows), sep = ""))
       print(paste("N = ", sum(temp$N_rows_NotEqual_mode),
                   " RecordIDs with mean N_rows = ",
                   mean((temp %>% filter(N_rows!=mode.N_rows))$N_rows),
                   "\nPrinting problematic RecordIDs now. Problems stored in `problems`",
                   sep = ""))
       #print(paste("N = ", sum(temp$N_Equal_mode[1]), " samples with nrow = ", temp$N, sep = ""))
       
       problems <- temp %>%
         filter(N_rows != mode.N_rows) %>%
         select(RecordID, N_rows, mode.N_rows) %>%
         left_join(meta.MEGA.T21_visit1_Celiac02, by = "RecordID") %>%
         select(FamilyID, RecordID, MEGA.IID, MEGA.LabID, N_rows, mode.N_rows) %>%
         unique()
       
       print(problems)
     }

[1] "Great, now all RecordIDs have the same number of rows of data."

If all of the following equal the same number of rows, then this dataframe now represents genotyping only once per individual

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID, MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(RecordID, MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(MEGA.LabID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(MEGA.IID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(RecordID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID, MEGA.IID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID) %>%
  unique() %>%
  nrow()

[1] 210

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(RecordID, MEGA.IID) %>%
  unique() %>%
  nrow()

[1] 210

print("If all of the numbers printed above equal the same number of rows, then this dataframe now represents genotyping only once per individual.")

[1] "If all of the numbers printed above equal the same number of rows, then this dataframe now represents genotyping only once per individual."

print("Good.")

[1] "Good."

Final check that all RecordIDs now have only 1 MEGA.IID after using the EXCLUDE_from_analysis column to filter

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(RecordID, MEGA.IID) %>%
  unique() %>%
  group_by(RecordID) %>%
  summarise(N = n()) %>%
  arrange(desc(N)) %>%
  left_join(meta.MEGA.T21_visit1_Celiac02, by = "RecordID") %>%
  filter(N > 1)


print("Good. All RecordIDs have only 1 MEGA.IID after applying the EXCLUDE_from_analysis filter.")

[1] "Good. All RecordIDs have only 1 MEGA.IID after applying the EXCLUDE_from_analysis filter."

Check for twins and other relateds

If there are relateds in the dataset after applying the EXCLUDE filter, first prioritize keeping Celiac cases, then randomly choose one IID per FamilyID to retain for analysis

meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID) %>% #, MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  group_by(FamilyID) %>%
  summarise(N_RecordIDs_per_FamilyID = n()) %>%
  ungroup() %>%
  arrange(desc(N_RecordIDs_per_FamilyID)) %>%
  filter(N_RecordIDs_per_FamilyID > 1)

relateds <- meta.MEGA.T21_visit1_Celiac02 %>% filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID) %>% #, MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  group_by(FamilyID) %>%
  summarise(N_RecordIDs_per_FamilyID = n()) %>%
  ungroup() %>%
  arrange(desc(N_RecordIDs_per_FamilyID)) %>%
  filter(N_RecordIDs_per_FamilyID > 1) %>%
  left_join(meta.MEGA.T21_visit1_Celiac02, by = "FamilyID") %>%
  select(FamilyID, N_RecordIDs_per_FamilyID, RecordID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, Celiac, Sex, Age_at_visit1) %>%
  unique() %>%
  #left_join(meta_comorb.022021, by = c("FamilyID", "RecordID", "MEGA.LabID"="LabID")) %>%
  group_by(FamilyID) %>%
  arrange(desc(Celiac)) %>%
  mutate(index_within_FamilyID = row_number()) %>%
  ungroup() %>%
  select(FamilyID, N_RecordIDs_per_FamilyID, index_within_FamilyID, Celiac, Sex, Age_at_visit1, everything()) %>%
  arrange(FamilyID, index_within_FamilyID)

relateds

relateds %>%
  select(FamilyID, index_within_FamilyID, EXCLUDE_from_analysis, Celiac, Sex) %>%
  unique() %>%
  group_by(FamilyID, Sex) %>%
  summarise(N = n())

`summarise()` has grouped output by 'FamilyID'. You can override using the `.groups` argument.

print("Decision: From FamilyID F0348 (2 individuals), keep INVTH328HA9 because this person is a Celiac case.")

[1] "Decision: From FamilyID F0348 (2 individuals), keep INVTH328HA9 because this person is a Celiac case."

print("Note: All other familys have only males or only females, except for one family with 1 of each.")

[1] "Note: All other familys have only males or only females, except for one family with 1 of each."

print("Decision: Therefore, from all other FamilyIDs, randomly select one individual to remain in the analysis.")

[1] "Decision: Therefore, from all other FamilyIDs, randomly select one individual to remain in the analysis."

set.seed(1234) # Set seed for randomly shuffling rows
relateds_to_keep = relateds[sample(1:nrow(relateds)), ] %>%
  group_by(FamilyID) %>%
  mutate(random_index = row_number()) %>%
  ungroup() %>%
  mutate(Keep = ifelse(FamilyID == "F0348" & Celiac == 1, 1,
                       ifelse(FamilyID == "F0348" & Celiac == 0, 0,
                              ifelse(random_index == 1, 1, 0)))) %>%
  select(FamilyID, RecordID, MEGA.IID, MEGA.LabID, N_RecordIDs_per_FamilyID, index_within_FamilyID, random_index, Keep, everything()) %>%
  arrange(FamilyID, index_within_FamilyID) %>%
  filter(Keep == 1)
relateds_to_keep

relateds_to_keep$MEGA.IID

[1] "WG1-DNA_F10_HTP0295A"  "WG1-DNA_H05_HTP0262A"  "WG3-DNA_C03_HTP0333A3" "WG2-DNA_G10_HTP0480A"  "WG3-DNA_C05_HTP0432A"  "WG3-DNA_F05_HTP0566A" 

Update the EXCLUDE_from_analysis and EXCLUDE_reason columns to indicate which IIDs to exclude within multiple-IID FamilyIDs

meta.MEGA.T21_visit1_Celiac03 <- meta.MEGA.T21_visit1_Celiac02 %>%
  select(FamilyID, RecordID) %>% #, MEGA.IID, MEGA.LabID) %>%
  unique() %>%
  group_by(FamilyID) %>%
  mutate(N_RecordIDs_per_FamilyID = n()) %>%
  ungroup() %>%
  arrange(desc(N_RecordIDs_per_FamilyID)) %>%
  mutate(Multi_IID_FamilyID = ifelse(N_RecordIDs_per_FamilyID > 1, 1, 0)) %>%
  left_join(meta.MEGA.T21_visit1_Celiac02, by = c("FamilyID", "RecordID")) %>%
  mutate(Keep_from_Multi_IID_FamilyID = ifelse(MEGA.IID %in% relateds_to_keep$MEGA.IID, 1, 0)) %>%
  mutate(EXCLUDE_reason = gsub("[.]", "", EXCLUDE_reason)) %>%
  #filter(Multi_IID_FamilyID==1) %>% select(FamilyID, RecordID, N_RecordIDs_per_FamilyID, MEGA.IID, Keep_from_Multi_IID_FamilyID) %>% unique() %>% arrange(FamilyID, desc(Keep_from_Multi_IID_FamilyID))
  mutate(EXCLUDE_from_analysis = ifelse(Multi_IID_FamilyID == 1 & Keep_from_Multi_IID_FamilyID == 0, 1, EXCLUDE_from_analysis),
         EXCLUDE_reason = ifelse(is.na(EXCLUDE_reason) & Multi_IID_FamilyID == 1 & EXCLUDE_from_analysis == 1, "FamilyID with multiple RecordIDs",
                                 ifelse(!is.na(EXCLUDE_reason) & 
                                          Multi_IID_FamilyID == 1 &
                                          EXCLUDE_from_analysis == 1,
                                        paste(EXCLUDE_reason, "; FamilyID with multiple RecordIDs", sep = ""), EXCLUDE_reason))) %>%
  select(-c(Keep_from_Multi_IID_FamilyID)) %>%
  rename(N_T21_RecordIDs_in_FamilyID = N_RecordIDs_per_FamilyID) %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Multi_IID_FamilyID, N_T21_RecordIDs_in_FamilyID, everything()) #%>%
  #mutate(EXCLUDE_from_analysis = ifelse(is.na(EXCLUDE_from_analysis), 0, EXCLUDE_from_analysis))

meta.MEGA.T21_visit1_Celiac03

meta.MEGA.T21_visit1_Celiac03 %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID,
         Multi_IID_FamilyID, N_T21_RecordIDs_in_FamilyID, Celiac, EXCLUDE_from_analysis, EXCLUDE_reason) %>%
  unique() %>%
  arrange(desc(Celiac), desc(EXCLUDE_from_analysis))

Verify that the EXCLUDE variable is working correctly

Not excluding any Celiac cases unless they represent a RecordID genotyped twice on MEGA. Producing a dataset with only one MEGA.IID per RecordID. Producing a dataset with only one RecordID per FamilyID.

meta.MEGA.T21_visit1_Celiac03 %>%
  select(MEGA.IID, MEGA.LabID, Celiac, EXCLUDE_from_analysis, EXCLUDE_reason) %>%
  unique() %>%
  group_by(Celiac, EXCLUDE_from_analysis, EXCLUDE_reason) %>%
  summarise(N = n())

`summarise()` has grouped output by 'Celiac', 'EXCLUDE_from_analysis'. You can override using the `.groups` argument.

print("Good, not excluding any Celiac cases unless they represent a RecordID genotyped twice on MEGA.")

[1] "Good, not excluding any Celiac cases unless they represent a RecordID genotyped twice on MEGA."

print("Reasons we are currently excluding IIDs represented in the plink *.fam are:
      1) If the IID is for a RecordID that got genotyped twice, or
      2) If the IID is from a FamilyID with multiple genotyped individuals.")

[1] "Reasons we are currently excluding IIDs represented in the plink *.fam are:\n      1) If the IID is for a RecordID that got genotyped twice, or\n      2) If the IID is from a FamilyID with multiple genotyped individuals."

meta.MEGA.T21_visit1_Celiac03 %>%
  filter(EXCLUDE_from_analysis == 0) %>%
  select(RecordID, MEGA.IID) %>%
  unique() %>%
  group_by(RecordID) %>%
  summarise(N = n()) %>%
  arrange(desc(N))

print("Good, only one MEGA.IID per RecordID.")

[1] "Good, only one MEGA.IID per RecordID."

meta.MEGA.T21_visit1_Celiac03 %>%
  #filter(EXCLUDE_from_analysis == 0) %>%
  select(FamilyID, RecordID, MEGA.IID) %>%
  unique() %>%
  group_by(FamilyID, RecordID, MEGA.IID) %>%
  summarise(N = n()) %>%
  arrange(desc(N))

`summarise()` has grouped output by 'FamilyID', 'RecordID'. You can override using the `.groups` argument.

print("Good, only one RecordID per FamilyID.")

[1] "Good, only one RecordID per FamilyID."

View who is being excluded at this stage and why

meta.MEGA.T21_visit1_Celiac03 %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, Celiac, EXCLUDE_from_analysis, EXCLUDE_reason) %>%
  unique() %>%
  group_by(EXCLUDE_from_analysis, EXCLUDE_reason, Celiac) %>%
  summarise(N = n())

`summarise()` has grouped output by 'EXCLUDE_from_analysis', 'EXCLUDE_reason'. You can override using the `.groups` argument.

meta.MEGA.T21_visit1_Celiac03 %>%
  select(FamilyID, RecordID, MEGA.FID, MEGA.IID, MEGA.LabID, EXCLUDE_from_analysis, EXCLUDE_reason, Karyotype, Celiac, Sex) %>%
  unique() %>%
  group_by(EXCLUDE_from_analysis, EXCLUDE_reason, Karyotype, Celiac, Sex) %>%
  summarise(N = n()) %>%
  ungroup() %>%
  #filter(#EXCLUDE_from_analysis == 0 &
  #         !is.na(Celiac)) %>%
  mutate(Total = sum(N),
         Pct = round(100*N/Total, digits = 2)) %>%
  mutate(`N (%)` = paste(N, " (", Pct, ")", sep = "")) %>%
  select(EXCLUDE_from_analysis, EXCLUDE_reason, Karyotype, Celiac, Sex, `N (%)`) %>%
  mutate(Group = ifelse(Celiac == 0, "T21 without CD",
                        ifelse(Celiac == 1, "T21 with CD", NA))) %>%
  select(EXCLUDE_from_analysis, EXCLUDE_reason, Group, Sex, `N (%)`)

`summarise()` has grouped output by 'EXCLUDE_from_analysis', 'EXCLUDE_reason', 'Karyotype', 'Celiac'. You can override using the `.groups` argument.

Save to files

setwd(dir)

Warning: The working directory was changed to /Users/shawjes/Dropbox/EspinosaGroup/ANALYSIS/Celiac_MultiOmics/GRS/DSMIG_Shared/Manuscript_Figure1/Data inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

fwrite(meta.MEGA.T21_visit1_Celiac03, "MEGA_041822_META_CeliacGRS_v0.1_JRS.csv")
fwrite(meta.MEGA.T21_visit1_Celiac03, "MEGA_041822_META_CeliacGRS_v0.1_JRS.tsv", sep = "\t")